Here's an example of the utility in action for a quarter degree grid:

Codecov Report

Attention: 119 lines in your changes are missing coverage. Please review.

Comparison is base (c760cf7) 2.03% compared to head (badf930) 7.81%.

@@           Coverage Diff            @@
##            main     #35      +/-   ##
========================================
+ Coverage   2.03%   7.81%   +5.78%     
========================================
  Files          9      12       +3     
  Lines        591     678      +87     
========================================
+ Hits          12      53      +41     
- Misses       579     625      +46     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

50th degree bathymetry for the mediterranean:

You'll need to remove connected regions, otherwise, when prescribing fluxes it might be that temperature (or salinity) grow beyond physical values because there is no way to mix isolated cells.

I found a very simple way to do it using the scikit package in pyhton

using PyCall

const ABOVE_SEA_LEVEL = 0.01
const scikitimage = PyNULL()
copy!(scikitimage, pyimport_conda("skimage.measure", "scikit-image"))


function remove_connected_regions(bat)

    bathymetry = deepcopy(bat)
    batneg     = deepcopy(bathymetry)

    batneg[batneg.>0] .= 0
    batneg[batneg.<0] .= 1

    labels = scikitimage.label(batneg, connectivity = 1)
    total_elements = zeros(maximum(labels))

    for i in 1:lastindex(total_elements)
        total_elements[i] = sum(labels[labels.==i])
    end

    ocean_idx = findfirst(x -> x == maximum(x), total_elements)
    second_maximum = maximum(filter((x) -> x != total_elements[ocean_idx], total_elements))

    # the bering sea is disconnected from the pacific/atlantic ocean if we go below 80 degree latitude
    bering_sea_idx = findfirst(x -> x == second_maximum, total_elements)

    labels = Float64.(labels)
    labels[labels.==0] .= NaN

    for i in 1:length(total_elements)
        if (i != ocean_idx) && (i != bering_sea_idx)
            labels[labels.==i] .= NaN
        end
    end

    bathymetry .+= labels
    bathymetry[isnan.(bathymetry)] .= ABOVE_SEA_LEVEL

    return bathymetry
end

Thanks @simone-silvestri . In the original post I called his

Eventually I would also like to add the kwarg fill_inland_seas.

By "filling inland seas", I mean the same thing you mean "remove connected regions" (I think). An inland sea is a region that isn't connected to the ocean.

As for the solution, I think we want to avoid depending on python packages in ClimaOcean if we can --- is there any other way?

Hmmm, I haven't found any simple way like this in Julia, but there might be an equivalent package to scikit-image

Sorry @glwagner, this slipped my head long now!

I heuristically found that smoother bathymetry allows larger time steps and a lower bathymetry noise.
It might be that a conservative regridding leads to a very unsmooth and jagged bathymetry with spurious grid-scale features.
In my personal experience, I found out that a Gaussian interpolation (akin to repeated bilinear interpolations in passes on progressively coarser grids) is great for retaining the large-scale features that can be represented while eliminating the small-scale features that just show up on the grid as grid-scale noise

@simone-silvestri have you done many simulations at 1 and 2 degree? I think the averaging becomes more important as coarse resolutions (ie to preserve basic properties like ocean volume and connectivity).

I don't think we need to have "only one" way to generate bathymetry either. But it's good to know that we may want other different utilities.

I do think we have to interpolate in general (because we don't have a general regridder, and we will often want to do simulations on non-latitude-longitude grids). So I was thinking that we could consider a two step generation: first conservatively map from 1/60th to a resolution a bit closer to what we are targeting. The second step is to interpolate the remapped bathymetry to a finer grid.